
TECHNICAL FIELD OF THE INVENTION 

This invention relates in general to methods and 
apparatus for determining the range from one aircraft to 
another aircraft, and more particularly relates to methods 
and apparatus for determining the range through use of 
passive energy detection with directional receivers. 
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BACKGROUND OF THE INVENTION 

The goal of air-to-air passive ranging is to 
determine the range (distance) from one aircraft, called 
"ownship", to another, called "target", by means of 
electromagnetic energy emanating from target. For this 
5 purpose, ownship is equipped with directional receivers 

operable to determine the angle of arrival of the incoming 
energy at various times during a data collection interval. 
In most applications, the energy detected will be in the 
radio frequency (RF) or infrared ( IR) portions of the 
10 spectrum. By combining the directional measurements with 

measurements of ownship' s position and heading, the range 
between the two aircraft may be determined. The dominant 
problem in air-to-air passive ranging is ill-conditioning. 
Ill-conditioning refers to small errors in the directional 
15 measurements causing much larger errors in the range 

estimate. The sources of ill-conditioning include the 
limited baseline over which data is measured and the need 
to infer target motion from the data. 

Range is estimated from bearing measurements at 
20 different locations using some variation of a 

"triangulation principle". A "baseline" is the distance 
between two points at which the bearing measurements, 
called sightings, are taken. Other factors being equal, 
sightings projected from a long baseline intersect at a 
greater angle than those taken from a short baseline. 
Since sightings projected from a short baseline are nearly 
parallel, a small error in the measured angle results in a 
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large error in determination of range. Hence, range 
measurements determined using a long baseline are much 
more impervious to noise contained in the bearing 
measurements than are those determined over a short 
baseline. However, the desire to determine range as 
quickly as possible restricts the baseline that may be 
developed during data collection. Furthermore, flight 
path geometry also dictates the maximum baseline which may 
be used; for example, if two aircraft are flying towards 
one another, no baseline will be generated regardless of 
time . 

In air-to-ground ranging, the target can be 
presumed as being fixed in space. Thus, it can be 
described by a static model consisting of its coordinates 
in 3-dimensional space. However, this is not true of the 
air-to-air problem, where target motion must be 
considered. The additional degrees-of-f reedom implied by 
the dynamic model must be resolved from the data. This 
aggravates ill-conditioning because it diminishes the 
redundancy by which information may be distinguished from 
noise in the data. Unresolved noise in the data produces 
ranging errors. 

As established by the foregoing, a need has 
arisen to provide a method and apparatus for measuring the 
range from one aircraft to another, capable of supplying 
ranging information quickly, while minimizing the effect 
of small errors in data measurement on the ranging 
determination. 
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SUMMARY OF THE IbA^ENTION 

In accordance with the present invention, a 
air-to-air passive ranging apparatus is provided which 
substantially eliminates or prevents the disadvantages and 
problems associated with prior art passive ranging 
5 systems. 

In accordance with another aspect of the 
^ invention, a ranging system/? is provided wherein model 

data estimating the flight characteristics of the target 
plane are produced and compared with actual flight path 
data of the target plane. An error measurement between 
the model data and the actual flight path data is 
calculated and adjustment is made to model data which 
minimizes the error measurement. 

In yet another aspect of the invention, the 
error measurement comprises a means squared sighting 
error, a velocity related error, and three position 
related errors. A system error is defined as the sum of 
the means squared sighting error, the velocity error and 
the three position errors. A perturbation model is 
determined such that the sum of the model data and the 
perturbation model data will result in a new set of model 
data which reduces the error measurement. 

In a further aspect of the invention, a method 
and apparatus for initializing the model data is provided. 
In yet a further aspect of the invention, a method and 
apparatus for determining a flight path which minimizes 
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the effect of noise on the sighting measurements is 
provided. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

For a more complete understanding of the present 
invention and the advantages thereof, reference is now 
made to the following descriptions taken in conjunction 
with the accompanying drawings in which: 

FIGURE 1 is a perspective view of an aerial 
encounter between two aircraft; 

FIGURE 2a-b is a diagram illustrating the 
ranging problem; 

FIGURE 3 is a block diagram of the ranging 
system functional architecture; 

FIGURE 4 is a flow chart describing calculation 
of the perturbation model; 

FIGURE 5 is a flow chart describing calculation 
of the perturbation magnitude; 

FIGURE 6 is a flow chart describing calculation 
of the initial brackets; 

FIGURE 7 is a flow chart describing a flight 
path advisor; 

FIGURE 8 is a diagram describing the formation 
20 of the initialization model; 

FIGURE 9 is a graph of the system eigenvalues; 

and 

FIGURE 10 is a circuit diagram of the ranging 

system. 
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DETAILED DESCRIPTION OF THE INVENTION 

The preferred embodiment of the present 
invention is best understood by referring to FIGURES 1-10 
of the drawings, like numerals being used for like and 
corresponding parts of the various drawings. 

FIGURE 1 illustrates an aerial encounter between 
two aircrafts. Ownship 10 desires to determine the range 
between itself and a second aircraft, target 12. Target 
12 emits electromagnetic energy, generally indicated by 
the reference numeral 14. Ownship 10 contains directional 
receivers 16 operable to determine the angle-of -arrival of 
the incoming energy forces 14 at various times during a 
data collection and interval. In addition to the bearing 
measurements, fixes of ownship' s position and heading are 
monitored by the on-board inertial navigation system (INS) 
15 18. A rectangular coordinate system 19 having an axis 

20, an axis 22, and an axis 24, has its center 26 

located at ownship. The coordinate system 19 moves in 
pure translation with ownship 10; however, the directions 
of the coordinate axes 20-24 are unrelated to the body 
axes of ownship 10 but are defined in inertial space by 
the INS 18. For convenience, the bearing measurements 
from the directional receivers 16 are converted to 
measurements related to the fixed axes 20-24. The 
necessary conversions are readily computed from knowledge 
of ownship 's heading as provided by the INS 18. 

The location of the target 12 can thus be 
described as a point in the coordinate system 19 expressed 
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in terms of the three axes 20, 22, and 24. Alternatively, 
the target can be described in terms of angles from two of 
the axes, denoted E (for elevation) and A (for azimuth) 
and p (range). It should be noted that an unconventional 
definition of elevation has been adopted as being the 
angle between the axis 24 and the vector from the 

center 26 to the target 12. The angle E and A can be 
computed from the bearing measurements provided by the 
directional receivers. 

The design of the directional receiver 16 
depends upon the spectrum of the energy 14 being 
monitored. For most applications, the energy will be in 
either the RF (radio frequency) or IR (infrared portions) 
of the energy spectrum. RF systems are designed to 
monitor radar transmissions and RF communications from the 
target 12. An RF system is capable of detecting energy 
from very long distances, substantially exceeding the 
distances at which the target could acquire ownship on its 
own radar. Since ownship 10 is using a passive ranging 
system, i.e., not transmitting energy signals itself, 
ownship 10 may be able to locate target 12 without 
reciprocal detection. This results in a primary advantage 
of passive ranging, namely, stealth. Furthermore, long 
acquisition range also affords ownship more time to 
respond to the target's presence. In a RF system, the 
directional receiver 16 may be implemented by using a 
2-axis RF interferometer, or multiple interferometer units 
covering different sectors. 
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A directional receiver 16 capable of detecting 
infrared energy signals 14 may be designed to respond to 
thermally emitted IR radiation from the target. Because 
of the shorter wave length IR energy signals, the 
direction of the incoming IR energy can be resolved more 
precisely. Consequently, the range between the aircraft 
can be determined more quickly. The infrared system has a 
further advantage in that performance is not linked to 
target transmission protocol. However, the maximum 
acquisition range of the IR system is much shorter than 
the RF system, thus allowing the target 12 the advantage 
of detecting ownship 10 on its radar or with a comparable 
IR system. A directional receiver 16 for the IR system 
may be implemented using a FLIR (Forward Looking Infra 
15 Red) imager and video tracker. 

The angles E and A are subject to measurement by 
the directional receiver 16; however, the range must be 
determined from the bearing measurements and the INS data. 
FIGURES 2a-b illustrate that the range between two 
20 aircraft 10 and 12 cannot be determined solely on the 

direction of the sighting vector 25, i.e., the angles E 
and A, in the absence of acceleration of ownship 10. 
FIGURE 2 shows the flight path of ownship 10. The 
positions of ownship are shown at t^, tj^, x.^ and t 
25 Vector 28 depicts the distance and direction of the flight 

path of ownship 10 between time periods t^ and t ; 
likewise, vector 30 indicates ownship flight path between 
t^ and t^, while vector 32 indicates ownship' s flight path 
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between time periods and t3 . Sight line 34 indicates 

the sight vector between ownship 10 and target 12 at time 
t^; similarly, sight line 36 indicates the sight vector at 
time t^, sight line 38 indicates the sight vector at time 
t2 and sight line 40 indicates the sight vector at time 
t^. Starting point 41 indicates the position of target 12 
at the time the first measurement is taken. Vectors 42, 
44 and 46 indicate the direction and distance traveled by 
target 12 during the time period from t^ to t3 . From 
phantom starting point 47, phantom target vectors 48, 50 
and 52 indicate a flight path which would provide a 
possible solution to the range algorithms. From another 
phantom starting point 53, phantom target vectors 54, 56, 
and 58 are vectors which would also provide a solution to 
the ranging algorithm. 

At time t^, ownship 10 would receive an energy 
signal from target 12. From this signal, the directional 
receivers 16 aboard ownship 10 would be able to determine 
the direction of the sighting vector 25 between ownship 10 
and target 12, but would not be able to determine the 
magnitude (range) of the sighting vector 25. Thus, at 
time tQ, ownship 10 can only identify target 12 as being 
somewhere along sighting line 34. From t^ to t^, ownship 
flies along vector 28. At time t^, ownship takes a 
another bearing measurement, rendering sight line 36. 
Once again, ownship 10 knows that target 12 is somewhere 
along sight line 36, but the exact whereabouts of target 
12 is yet unknown. From time t^ to time t^, ownship flies 
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along vector 30. Since ownship is flying at a constant 
velocity and direction, vectors 28 and 30 are equal in - 
magnitude and colinear. At time t^ , another bearing 
measurement by ownship 10, resulting sight line 38. 
Whereas the range still cannot be determined from the 
sight vectors, some pertinent information can be 
determined from the data. 

^irst., it must be assumed that target 12 is 
flying at a constant velocity and in a constant direction. 
This is a valid assumption since most aircraft fly at a 
nominal speed. Given the assumption of constant velocity, 
the distance travelled by target 12 between t^-t^^ and 
^1"''^2 '"^s^ equal. Also, assuming constant direction, 

the vectors representing the flight path of target 12 must 
15 colinear. Given the sight lines 34-40, as shown in FIGURE 

2a, only one set of vectors will meet the criteria of 
colinearity and equal magnitude for any given point along 
the sight line 34 at tg . For example, given the starting 
point 41 of target 12 along sight line 34, only the flight 
path given by vectors 42 and 44 will provide a solution in 
which the vectors are colinear and equal in magnitude. 
Likewise, given a starting point 47, only the flight path 
represented by vectors 48 and 50 will result in a solution 
having colinear vectors of equal magnitude. Further, 
25 given starting point 53 on sight line 34, only the vectors 

54 and 56 will provide a suitable solution. 
Unfortunately, in an actual sighting, the starting point 
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41 of target 12 will be unknown and must be determined by 
the data. 

Even after another time period, the ambiguity 
cannot be resolved, as illustrated by the sight path 
vector 46 of the actual target and vectors 52 and 58 of 
the phantom target. In order to resolve this ambiguity, 
ownship must accelerate or decelerate, as shown in FIGURE 
2b. 

FIGURE 2b depicts the same information as FIGURE 
2a, with the exception that ownship 10 decelerated after 
time period resulting in a reduced distance travelled 
over vector 32. Consequently, the direction of sight line 
40 taken at time t3 has changed from FIGURE 2a. Due to 
the deceleration of ownship 10 at X.^. without a 
corresponding deceleration by target 12, only the flight 
path defined by vectors 42, 44 and 46 meet the criteria of 
colinearity and equal magnitude, and hence, the ambiguity 
is resolved. Thus, the magnitude of the sight vector can 
be determined from the angles of the sight lines 34, 36, 
20 38 and 40. 

However, as can be seen from FIGURE 2a, the 
correct determination of the range is highly dependent 
upon a accurate measurement of the sight line. For 
example, a slight change in the angle of sight line 40 
25 taken at t3 in FIGURE 2b could result in the intersection 

of the sight line 40 with the endpoint of vector 52, 
rather than the endpoint of the actual target vector 46. 
As can be seen in FIGURE 2b, this small error in angle 
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measurement could result in a miscalculation in range of 
greater than 20%. 

The effects of ill-conditioning may be reduced 
by introducing data redundancy into the equations used to 
determine the range, i.e., solving the range algorithm 
over several additional sightings. However, data 
redundancy may be counterproductive if achieved by 
reducing the sample interval. 

The inability of the direct solution method to 
deal with ill-conditioning has been overcome by the 
present invention using a method which does not compute 
range directly, but rather computes the parameters of a 
kinematic model of the target. Once an acceptable model 
is created, range may be computed from the model, rather 
15 than the actual bearing measurements. Since the model has 

a fixed number of parameters, redundancy increases 
monotonically as more data are assimilated. 

FIGURE 3 illustrates a functional organization 
of the air-to-air ranging system functional architecture. 
Input data 60 is stored in a measurement buffer 62 which 
holds the input data 60 for successive time frames. Data 
from the measurement buffer 62 is used to create the 
initialization model 64. Initialization model 64 is 
available to a buffer containing the resident model 66 
25 from which range estimates 68 are computed. Resident 

model 66 is also used to generate estimated azimuth and 
elevation "synthetic measurements" 70. The synthetic 
measurements 70 are compared with data from the 
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measurement buffer 62 in a subtracter 72 yielding sighting 
errors 74. The sighting errors 74 are used in conjunction 
with the resident model 66 to construct sighting error 
model 76. A position penalty model 78 is determined using 
data from the measurement buffer 62, the resident model 
66, and data regarding the maximum acquisition range 80. 
A speed penalty model 82 is assessed using the resident 
model 66 and data regarding the nominal target speed 84. 
The position penalty model 78 and the speed penalty model 
82 are added to the sighting error model 76 in summation 
block 86. Data from the summation block 86 along with 
data from the position penalties 78 and the speed penalty 
82 are used to construct an system error model 88. The 
system error model 88 is used to produce a perturbation 
15 model 90. The perturbation model 90 updates the resident 

model 66. Since a single perturbation of the resident 
model 66 may not provide a sufficient accuracy between the 
resident model 66 and the input data 60, the algorithm is 
iterated until the resident model 66 describes the input 
20 data 60 within a predetermined range. After the resident 

model 66 is sufficiently accurate to correspond with the 
input data 60, new input data 60 stored in the measurement 
buffer 62 relating to the next time period is compared 
with the data from the resident model 66 and the process 
25 repeats itself. 

In operation, the input data 60 is received from 
the IMS 18 and directional receiver 16. The data includes 
measurements of time (t), azimuth (A), elevation (E), 
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ownship position relative to initial starting point (a^, 
^3), and ownship heading (y^, y^ and y^ ) . Azimuth and 
elevation data are produced by the directional receivers 
16 and are adjusted to conform to the coordinate system 
19. Ownship position and heading are produced by the INS 
18 and likewise conform to the coordinate system 19. 
Timing data is received from an internal clock. 

The input data 60 is stored in a measurement 
buffer 62. While the input data can be processed 
substantially in real-time, the measurement buffer 60 is 
necessary to store the data for use in different tasks, 
such as producing the initialization model 64, the 
sighting error 74 and the position penalty 78. 

The initialization model block 64 computes the 
initial estimate of the model parameters. The six model 
parameters, denoted as the vector m, describe the position 
of the target 12 as a function of time, including three 
coordinates of position (x^^, vi^ and X3 ) and a three 
dimensional velocity vector having components v^^, Vj and 
20 ^3- Since target 12 is assumed to be flying at a constant 

speed and direction, the velocity y can be assumed 
constant. Calculation of the initialization model block 
64 will be described in greater detail in conjunction with 
FIGURES 8 and 9. 

25 Once the initial model 64 has been determined, 

it is used as the first resident model 66. The resident 
model 66 is the system's current best estimate of the 
target location and velocity. The ranging algorithm 
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revises the resident model 66 at each time interval to 
correspond with the actual input data over a period of 
time. Once the resident model 66 can no longer be 
improved, the system's best estimate of range can be 
obtained from the model parameters. As each new set of 
input data 60 is received, a perturbation model 90 is 
formed to update the resident model 66. The computation 
of the perturbation model is discussed in connection with 
FIGURES 3-6. 

Synthetic measurements 70 are measurements of 
elevation (E) and azimuth (A) calculated from the data 
from the resident model 66. The synthetic measurements 70 
would correspond to the actual input data 60 if the 
resident model 66 was an accurate description of the 
15 target 12 and the input data 60 was accurate and 

noise-free. The relationship between the azimuth and 
elevation and the model parameters are as follows: 
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^m-i = tan 
mi 



X2(0) . t.v^ * a^Ct.) 

-1 



20 ^1<0) * ^i^i * ^i(ti) 
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^m-i = tan 
mi 



-1 



Xi(0)+t.v^ + a^(t. ))2 Mx2(0)*t.V2*a2(t. ))■ 



X3(0) . t.V3 . a3(t.) 
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where and are the values computed from the resident 

model 56 at time t. . 

1 

An error function J(m) is defined by the 



equation : 



N-1 



= 1/N Z (E^, - E„, )^ + (A . - A . )^ 



'oi -mi' <^oi " ''^mi) 

i=0 



where E^^ and A^^ are the observed elevation and azimuth 
10 at time t^ and E^^ and A^^ are the values computed from 

the resident model 66 ac time t^ . J(m) is the mean 
squared sighting error evaluated with the data from the 
resident model, m. The object of the ranging system is to 
determine a model which minimizes J(m) . Once this model 
15 is obtained, range is readily computed from the model 

parameters. To obtain the model which minimizes J(m) , the 
resident model 66 is "perturbed" in a direction and 
magnitude which produces a new resident model which 
improves correlation between the actual observed data 
stored in the measurement buffer 62 and the synthetic 
measurements 70 calculated from the model. 

The subtracter 72 provides the sighting errors 

^oi " ^mi ^oi ' Vi ' calculation of 

the mean squared error function. The mean squared 
25 sighting error along with the sighting error perturbation 

model 76. 

The mean squared sighting error of the perturbed 
model is denoted: 
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J' = J{m + 6m), where 6in is the perturbation. 
This equation (the "perturbation equation") may 
be approximated by the second degree model: 
J' = J(m) + Sm'^g + (1/2) 6m'^H 6m 

where 3 is the gradient of J(m), a vector whose 
elements are the partial derivatives of J(m) with respect 
to the model parameters, i.e.: 

9 J. = (aj/3m^), for r = 1, 2,... 6 

The vector 3 can also be described in relation 
to the sighting error as: 

N-1 

Sj = (-2/N) r (E.VE„, + a.7A .) 

1 mi 1 ma ' 

i=0 

where and are the instantaneous sighting errors for 
elevation and azimuth: 

~ ~ ^nd a. = A . - A . 

1 01 mi 1 01 mi 

and VE^ and VA are the gradients of E„ . and A . with 

* mi mi 

respect to the parameters of the model m. 

The "Hessian" matrix, H, is denoted as: 

25 H = [h^gl = ( (a^J)/Om^3mg) J 

where m^. and m^ respectively refer to the r^^ and s^^ 
parameter in the set of six model parameters. The H 



15 



20 



19 



10 



matrix for the sighting error may be approximated by the 
equation: 

N-1 

H = (2/N) I (VE^. VE^,'^ + VA . 7A 
J mi mi mi mi ' 

i=0 

Since the vectors involving VE . and 7A . have 
, , , mi 
already been computed in determining the gradient of J(m) 

(denoted as _g in the equations), no new differentiations 
are required to compute the H matrix. 

The Hessian matrix, H^, and the gradient of 
^(m)/ 3j/ of the sighting error perturbation equation 
J(m) are summed with corresponding Hessian matrices 
15 and gradients found in relation to speed and position 

penalties 78 and 82, described below, to render a system 
error model 88. The system error model 88 is used to 
create a perturbation model 90. 

Position penalties model 78 and speed penalty 
model 82 provide constraints to the ranging algorithm. 
"Constraints" provide a mechanism for including real-world 
knowledge which is not explicitly present in the input 
data 60 into the ranging algorithms. Constraints restrict 
the domain of feasible solutions to target models with 
25 realistic parameters. Four constraints are employed in 

the system: a speed constraint, and three position 
constraints. The position constraints include a maximum 
range constraint, an elevation position constraint, and an 



20 



20 



azimuth position constraint. The use of constraints 
speeds convergence of the resident model 65 and reduces 
the effect of ill conditioning. While the use of 
constraints requires additional computational costs, the 
5 costs are offset by reduction in the iterations required 

to achieve convergence of the resident model 66. 

The actions of the constraints are indirect 
rather than direct; since the constraints control 
functions of the model parameters rather than the 
10 parameters themselves. For instance, the speed constraint 

acts upon the magnitude of the target velocity, i.e., 
^^1 * ^2 * V3^)*5. Clearly, the speed may remain 
constant even though two or more of the velocity 
parameters change. Also, it is clear that changes in the 
15 position parameters, x^, y.^ and X3 do not effect the speed 

of the target. Therefore, all target parameters may 
change without affecting the speed, and hence, the speed 
constraint. 

The speed constraint imposes a penalty when the 
20 velocity parameters, v^ , v^ and V3 , describe a target 

which has a magnitude of velocity deviating significantly 
from a nominal speed. Most aircraft of the same class fly 
at approximately the same speed, thus a predetermined 
nominal speed may be set, and a speed constrained 
25 perturbation model may be determined that requires that 

the target be within a specified tolerance of the 
predetermined nominal speed. 
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The function used to determine the speed 
penalty is : 

^ = K[ (y.yj^)/Ay]^^ 

where: y = target speed determined from model 
Pj^ = nominal (predetermined) speed 

K = weight of the constraint 

Ay = bandwidth of the constraint 

M = order of the constraint 



All of the elements of the speed penalty 
function are non-negative numbers. The element "y" is 

15 usually determined as the magnitude of the velocity 

factor, while is chosen as the nominal speed of a 

particular aircraft. The bandwidth. Ay, is choosen in 
proportion to the uncertainty with which the true target 
speed is known. The weight = K, is chosen equal to the 

20 desired value of the speed penalty when model speed 

deviates from nominal speed by one bandwidth. The order, 
M, increases the influence of the constraint beyond the 
bandwidth. In the preferred embodiment, orders above M = 
8 are not required, and are avoided, because of the 
25 associated numerical overflow and underflow problems in 

the evaluation of the derivatives in the stop band. 

The first and second order partial derivatives 
of lif are computed as constituents of the perturbation 
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equations described in conjunction with the sighting error 
perturbation model J(m) to render a corresponding Hessian 
matrix and a gradient vector 3^. Since the constraint 
is entirely velocity- related, the arguments are reduced 
from the six used in the sighting error perturbation 
equations, to three used in the speed constrained 
perturbation equations. Thus, in computing in the Hessian 
matrix of the speed constrained perturbation equations, 
only six second order partial derivatives are computed, 
rather than the twenty-one computed in connection with the 
sighting error perturbation equation. Furthermore, the 
constraint ^ is not time-dependent, as is the calculation 
for mean squared sighting error J(m). Therefore, 
time-iterative calculations need not be performed to 
15 compute ip. 

The position constraints require that the 
components of the target model associated with the initial 
position lie within a "cone" centered about the initial 
sighting vector. The angular cross-section of the cone is 

20 sized according to the uncertainty of sighting 
measurements, and is a statistical choice. 

A penalty is associated with each position 
constraints: a elevation position penalty, an azimuth 
position penalty, and a maximum range penalty. The 

25 equations for these penalties are: 



^ = ^ f<Eoi-Emi>/^E) 
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Pj. = exp [a(p - 9^)/9^] 

The constants K^, K^, K^, Me, Ma, and a are 
chosen as described in conjunction with the speed penalty. 
AE and AA are chosen in proportion to the uncertainty with 
which the true target elevation and azimuth are known. 
The factor is set equal to the maximum acquisition 

range. 

Three position penalty models are constructed, 
one for each of the three penalty equations. The position 
penalty models are determined as described above in 
connection with the sighting error model rendering 
equations for a Hessian matrix and gradient vector for 
each penalty model. As with the velocity penalty, the 
Hessian matrix for the position penalty equations is 
limited to the three position parameters - x^, and x^ . 

After computation of the sighting error model 
76, speed constrained model 82 and three position models 
78, the models are summed to create a system model 88. 
The equations used to define the system model are: 

J_ = J(m) +*+P +p +p 
s ear 

9s - 3j ^ 3^ - * 9pa * 9pr 

Several methods exist for minimizing the mean 
squared error Jg(m). The preferred method is the 
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Newton-Raphson method, which converges faster and provides 

for larger perturbations than other methods, such as the 

gradient method. A key step of the Newton-Raphson method 

is to solve 6m = 3^, A flow chart of the preferred 

5 method for solving this equation is given in FIGURE 4. In 

block 92, a Cholesky factorization is performed. The 

Cholesky factorization involves factoring the H matrix 

s 

into lower and upper diagonal matrices, L and l"^. In 
block 94, the equation = 3^ is solved for the vector x 

10 (L and 3^ are known). In block 96, the equation, 6m = y 

is solved for 6m. 

The accuracy of the perturbation model, 6m, 
computed in this manner, may be improved as illustrated in 
blocks 98-107 of FIGURE 4. In block 98, the solution for 
6m from block 96 is set as the initial solution, 651^. In 
block 100, a residual, r^, is determined from the equation 
_ T 

2 - L L fiiDo' 'Jsing a double precision computer 
arithmatic. In block 102, the equation L L^Yq ' 
20 solved for the "correction", y^. In block 104, y^ is 

added to 6mQ to provide a more accurate solution for the 
perturbation, 6m. As indicated in decision block 106, if 
more accuracy is desired the loop starting at block 98 
may be repeated. This algorithm used to improved accuracy 
of the initial perturbation solution has the advantage 
requiring only one double precision calculation for the 
remainder in block 100; the other two calculations in 
blocks 102 and 104 may be performed using single precision 
computations. 
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The system error model has both direction and 

magnitude. However, to establish the optimum magnitude of 

the perturbation, the preferred method is to perform a 

linear search, which allows non-linearities of higher 

degree in the system error function to influence the 

optimum magnitude. The purpose of the linear search is to 

determine the weight of the pertubation 6 at which the 

— m 

system error is minimized. In the preferred embodiment 
Davidon's cubic interpolation method is utilized. 
Davidon's cubic interpolation method is described in 
detail in "Methods of Optimization", G.R. Walsh, John 
Wiley & Sons (1975). 

A flow chart outlining the Davidon algorithm is 
illustrated in FIGURE 5. In block 108, a "bracket" having 
an upper and lower boundaries is established. Within the 
bracket, the optimum weight (magnitude) of the 
perturbation is located. The bracket may be determined by 
a trial-and-trial procedure, or by the procedure discussed 
hereinafter in connection with FIGURE 6. 

Block 110 is the first step of iterative loop 
used to narrow the boundaries of the bracket. In block 
110, the width of the bracket, i.e., the distance between 
the upper and lower boundaries, is determined. In block 
112, the width of the bracket is compared to a 
25 predetermined value. The predetermined value corresponds 

to a width at which further precision in determining the 
optimum weight of the perturbation no longer warrants 
additional computational time. If the width of the 
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bracket is less than or equal to the predetermined value, 

then the optimum weight is set to a value equal to the 

arthimatic average of the boundary weights as illustrated 

in block 114. If the widths of the bracket is greater 

than the predetermined value in block 112, then Davidon's 

cubic interpolation algorithm is used to estimate the 

location of the optimum weight within the bracket in block 

118. Davidon's cubic interpolation method computes a 

cubic polynomial estimating the system error J based on 

s 

the values of the system error and the derivative of the 
system error at the upper and lower bracket boundaries. 
The minimum value of the cubic polynomial can be 
determined from its derivative. The estimated optimum 
weight is set to the point at which the cubic polynomial 
15 is at a minimum. In block 120, the actual system error 

and its derivatives in the direction of the perturbation 
are computed for the estimated optimum weight determined 
in block 118. In decision block 122, if the derivative 
determined in block 120 is less than zero, then the lower 
boundary of the bracket moves to the location of the 
estimated optimum weight determined in block 118, thus 
narrowing the width of the bracket. If the derivative is 
not less than zero in decision block 122, and the 
derivative is greater than zero in decision block 126, the 
upper boundary of the bracket moves to the location of the 
estimated optimum weight, as described in block 130. If 
the derivative is either less than zero in block 122 or 
greater than zero in block 126, then the loop repeated 
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starting at block 110. If the derivative is equal to zero, 
then the optimum weight is set equal to the estimated 
optimum weight. 

After a number of iterations, the width of the 
bracket will be less than the predetermined value given in 
block 110. At this point, the linear search is completed, 
and the optimum weight is selected as the average between 
the boundary weights in block 114, or as specified in 
block 128. 

In some instances, may be substantially flat 
^ around it5 minimum, resulting in a slow convergence. If 

successive iterations of the bracketing algorithm are 
producing new boundaries for one side of the bracket 
without alternation between upper and lower boundaries, 
15 then steps must be taken to aid convergence. For example, 

if successive iterations of the algorithm have produced 
only incremental changes in the upper boundary, then the 
lower boundary should be moved closer to the upper 
boundary by a predetermined factor, such as by one third 
of the width between the upper and lower boundaries. 

FIGURE 6 illustrates a flow chart for 
determining the initial bracket boundaries. In block 132, 
the system error and its derivative along the direction of 
the perturbation is computed for a weight equal to zero. 
In decision block 134. the derivative of the system error 
is checked for a negative condition; if the derivative is 
equal to or greater than zero, then an error flag is 
returned in block 136. indicating that no reduction of 
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system error is possible in the direction of the 
perturbation as calculated. If the derivative is less 
than zero, then the lower boundary is set for a weight of 
zero in block 138. Assuming a non-negative derivative of 
the system error, a first weight (W^) of the model 
perturbation is determined such that the larger of the 
velocity penalty or the position penalty will equal the 
system error determined in block 132. In block 142, the 
upper boundary is set to equal . In block 144, the 
derivative of the system error along the direction of the 
perturation at the upper boundary is evaluated for use in 
the Davidon's algorithm of FIGURE 5. I f the derivative of 
the system error is negative at , Davidon's cubic 
polynomial will not be able to compute a suitable 
approximation of the system error between the upper and 
lower boundaries. In this instance, the boundaries are 
relocated as described in block 145. If the value of the 
system error at exceeds or equals the value of the 

system error at the lower boundary, then the upper 
boundary is set at 1/3 . Otherwise, the lower boundary 
is set at W^^ and the upper boundary is set at 2Vi^. 

After the optimum weight of the perturbation is 
determined, the resident model 66 is updated by adding the 
perturbation multipled by the optimum weight to the 
resident model 66. The new resident model is tested to 
determine whether further iterations would be fruitful, 
i.e., whether successive perturbations of the resident 
model will have a significant effect on the range 
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determination. By calculating a magnitude of the 

perturbation, and comparing it with a predetermined limit, 

the capacity for further impro-vement of the resident model 

66 is evaluated. If no further iterations are necessary, 

the resident model is considered valid and can be used to 

compute the range at any time in the data collection 

interval. The system will continue iterate after the 

range determinations are valid in order to update the 

range information and to detect path deviations by taraet 
10 12. 

As described in connection with FIGURES 2a-b, it 
is necessary ownship 10 to accelerate relative to target 
12, in order to allow computation of the range. A method 
for determining the "optimum maneuver" which will produce 
15 acceleration (or deceleration) in a direction having 

maximum insensi ti vi ty to errors in the measured sighting 
vectors 25 is illustrated in the flow chart of FIGURE 7. 
In Block 146, the sighting vector (w) 25 is shown in terms 
of the elevation (E) and the azimuth (A). In Block 148, a 
matrix W is defined as the sum of the sighting vector 
times its transpose over N sightings. The eigenvalues of 
the matrix W are determined in Block 150. In block 152, 
an eigenvector, v, associated with the smallest eigenvalue 
found in Block 150, is determined. This eigenvector is 
used in Block 154 to advise ownship 10 of the change of 
ciirec-ion which will minimize errors due to noise in the 
sighting data. Hence, if ownship 10 is travelling with a 
normalized velocity vector of (1,0,0) and the 
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eigenvector, u, has a value of (0,1,0), the new direction 
would be equalled to v + ku = (l,k,0), where k is a 
positive or negative real value. The factor k is a 
predetermined value' chosen to reflect ownship's ability to 
5 change directions substantially instantaneously. The 

eigenvalues and eigenvectors of matrix W may be determined 
using well known mathimatical techniques. 

As previously shown, range cannot be determined 
when both target 12 and ownship 10 fly with constant 
10 velocity. However, range can be estimated by assuming 

target speed (the magnitude of the target's velocity). By 
determining the estimated range, all of the model 
parameters are established. Thus, the initialization 
model 64 may be established after an estimated range is 
15 computed. The preferred method of computing the 

initialization model 64 is illustrated in FIGURES 8 and 9 . 
The method uses the results obtained in the flight path 
advisor illustrated in FIGURE 7. 

FIGURE 8 illustrates relative target positions 
20 at time t^ and time t . . The axis u^ and are 

eigenvectors of the matrix W, defined in relation to 
FIGURE 7, which correspond to the largest and second 
largest eigenvalues respectively. Typical eigenvalues of 
the matrix W are illustrated in FIGURE 9. 
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To solve for the initial value of the target 
velocity, v_: 

= t<u^, Vq> - cos ((ej^/2) + ♦)] + 

where : 



10 ^1' ^2' -3 " eigenvectors corresponding to 

largest, second largest, and third largest 
eigenvalues of matrix W 
Vq - ownship velocity vector 

' angle between 9^ and 9^, calculated as 

15 H 

{12\^) ^, where is the second largest eigenvalue of 

matrix W. 

0 - angle between w^ and v^, where v is the 
relative velocity vector 
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The acquisition range is given by: 
Pq = [-B-(b2-AC)^]/A 



where A = 3 



B = Pf^u^, Vq> sin ((9j^/2) +0). 

_ _ 2 2 

C - " ^^0 ^T 

speeds of the ownship 

and target) 
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and 3 = ((f^ + + S^^jj.i 
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where = tan (0^); = tan ( 8^ ) 
f2 = [(t^/t^)^ (1 ^ s^2j2 _ 

2(t^/t2) (l^S^^) (S^2^(S^/S2)) 

the initial target position, c in FIGURE 9, can be 
determined as: 

c = Pq (cos( (8j^/2)u^ - sin (8j^/2)u2) 

Referring now to FIGURE 10, a hardware 
implementation of the ranging system is illustrated. A 
processor 180 is connected to a program memory 181 and a 
storage memory 182 and two arithmetic coprocessors 184 and 
186. Input data measurements 60 are input to the processor 
180, and range determination outputs are output from the 
processor 180. 

Comparing FIGURE 3 with the hardware 
implementation of FIGURE 10, the storage memory 182 is 
used to implement the measurement buffer 62 and to store 
interim calculations of the synthetic measurements 70, 
sighting errors 74, sighting error model 76, position 
penalty model 78, speed penalty model 82', system error 
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model 88, perturbation model 90, resident model 66, and 
initialization model 64. The processor 180 manages data 
transfers between the storage memory 182, the arithmetic 
coprocessors 184 and 186, and the data input and output 60 
and 188. The program memory 181 supplies the instruction 
for processor 180. The program memory 181 may be internal 
to the processor 180. In the preferred embodiment, a TMS 
320 signal processor is used as the processor 180, 
although other microprocessors may be substituted as well. 
The arithmetic coprocessors 184 and 186 are used to 
facilitate calculation of the Hessian matrices and the 
gradient vectors associated with sighting error model 76, 
the position penalty model 78, and the speed penalty model 
82. Also, the arithmetic coprocessors 184 and 186 are 
used to determine the magnitude and direction of the 
perturbation model and to calculate the initialization 
model 64. 

Input data measurements 60 are preferably 
filtered prior to input to the memory 182. The filter 
serves to enhance the signal to noise ratio of the data, 
while providing a down-sampling function. 

It is important to note that other 
implementations for the ranging system may be substituted 
for the implementation illustrated in FIGURE 10. For 
25 example, a plurality of processors 180 may be utilized in 

order to increase the processing speed of the 
implementation. Since increased processing speed may 
allow the input data to be received at shorter time 
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intervals, this alternative may be desirable for some 
applications. If multiple processors are used, it may be 
advantageous to provide separate memories for some or all 
of the processors in order to prevent a bottleneck from 
occurring at the storage memory 182. Alternatively, if 
greater processing speed is required, the ranging system 
could be implemented using a pipeline configuration, 
wherein discrete arithmetic components replaced the 
functions performed by the processor 180. However, while 
a pipeline configuration can greatly increase processing 
speed, the number of additional components used in this 
implementation may result in a system which is too large 
or too costly in comparison with the benefits of increased 
processing speed. 

■"■^ Although the present invention has been 

described in detail, it should be understood that various 
changes, substitution and alterations can be made therein 
without departing from the spirit and scope of the 
invention as defined by the appended claims. 
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TECHNICAL ADVAMTAGES OF THE I^?VENTION 

It is an important technical advantage of the 
invention that the range between two aircraft can be 
accurately calculated from data subject to noise. 

It is a further technical advantage that the 
invention can determine the range between two aircraft 
from passively detected data. 

It is another technical advantage that a flight 
path optimizing ranging performance may be calculated. 

It is yet another technical advantage that the 
invention supplies ranging information quickly. 

It is yet a further technical advantage that the 
invention provides range information from directional 
data. 

It is yet a further technical advantage that the 
15 invention provides a method for determining initial 

estimates of the targets position and velocity. 
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